A symmetry-breaking mechanism for the Z4 general-covariant evolution system 
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The general-covariant Z4 formalism is further analyzed. The gauge conditions are generalized with 
a view to Numerical Relativity applications and the conditions for obtaining strongly hyperbolic 
evolution systems are given both at the first and the second order levels. A symmetry-breaking 
mechanism is proposed that allows one, when applied in a partial way, to recover previously proposed 
strongly hyperbolic formalisms, like the BSSN and the Bona-Masso ones. When applied in its full 
form, the symmetry breaking mechanism allows one to recover the full five-parameter family of first 
order KST systems. Numerical codes based in the proposed formalisms are tested. A robust stability 
test is provided by evolving random noise data around Minkowski space-time. A strong field test is 
provided by the collapse of a periodic background of plane gravitational waves, as described by the 
Gowdy metric. 

PACS numbers: 



I. INTRODUCTION 

The waveform emitted in the inspiral and merger of a 
relativistic binary is a theoretical input crucial to the suc- 
cess of the laser interferometry gravitational laboratories 
P, 0, 0, 01 ■ Although the regular orbiting phase can be 
treated with good accuracy by well-known analytical per- 
turbation methods, the later phases belong clearly to the 
strong field regime of either a black hole or a neutron star 
collision, so that a computational approach is mandatory. 
This kind of computational effort has been the ob^iective 
of the Binary Black Hole (BBH) Grand Challenge and 
other world wide collaborations. The resulting numerical 
codes are based on the so called ADM formalism 6] for 
the Einstein field equations, where only a subset of the 
equations are actually used for evolution whereas the re- 
maining ones are considered as constraints to be iinposed 
on the initial data only (free evolution approach 

It is clear that, by taking the constraints out of the 
evolution system, one is extending the solution space. 
This extension is the crucial step that opened the way 
to new hyperbolic formalisms after the seminal work of 
Y. Choquet-Bruhat and T. Ruggeri j^, using the con- 
straints to modify the evolution s^tem in many different 
ways i Ea m d IllJliyaiia [H m, even taking 
additional derivatives [3, ■ These formalisms can 

be interpreted as providing many non-equivalent ways of 
extending the solution space of Einstein's equations with 
at least one common feature: constraint equations are 
left out of the final evolution system. 

This means that the resulting systems do have an 
extended solution space, which includes constraint- 
violating solutions in addition to the ones verifying the 
original Einstein's equations. As far as the constraints 
are first integrals of the extended evolution system, Ein- 
stein's solutions could be computed by solving the con- 
straints equations only for the initial data (free evolu- 
tion). But, unless some enforcing mechanism is used 
during the subsequent time evolution, numerical errors 



will activate constraint-violating modes. Numerical sim- 
ulations can deal with such modes, at least when the 
deviations from Einstein's solutions are moderate. But 
it happens that large deviations are usually associated 
with instabilities of constraint-violating modes, leading 
to code crashing. 

In particular, numerical codes based on these new hy- 
perbolic formalisms happen to be quite intolerant to vio- 
lations of the Hamiltonian constraint in the initial data. 
This is a serious drawback if one is planning to use the 
results of the analytical approximation to the regular or- 
biting phase of a binary system as initial data for a nu- 
merical simulation of the final ringdown and merger. The 
numerical code will crash before any template for the 
gravitational wave emission could be extracted. This is 
because the analytical data are just a good approxima- 
tion, so that the energy constraint does not hold exactly 
and the code is intolerant to that kind of "off-shell" initial 
data. Although there can be other options, we claim that 
a numerical code tolerant to constraint violation would 
be undoubtedly the best alternative. 

There have been some attempts in that direction. In 
Ref. "21] the technique of Lagrange multipliers was used 
for including the constraints into the dynamical system. 
The extended system includes the Lagrange multipliers 
as additional dynamical fields (A-system). A further step 
along that direction is given in Ref. 22], where the ex- 
tended system is "adjusted" by further combining the 
evolution equations with the constraints, including then 
many extra arbitrary parameters. In both cases, the goal 
is to enforce the constraints in a dynamical way. One can 
monitor the errors by looking at the " subsidiary system" , 
which can be derived from Bianchi identities and can be 
interpreted as the evolution system for constraint devi- 
ations. Parameters are adjusted in a way such that the 
characteristic speeds of the subsidiary system are either 
real and non-zero (so that the corresponding errors prop- 
agate away to the boundaries) or they have the right sign 
in the imaginary part to enforce damping (instead of ex- 
ploding) the constraint deviations. 
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A simpler option (but not the only one, see for instance 
[2^) for including the constraints into the evolution sys- 
tem would be the general-covariant extension of the Ein- 
stein field equations proposed recently (Z4 system) (25|: 



Rt^u + V,,Z, + V^Z,, = 8 TT (T^, - - T g^,) , (1) 

so that the full set of dynamical fields consists of the 
pair {g^^,Z^}. The solutions of the original field equa- 
tions can then be recovered by imposing the algebraic 
constraint 



(2) 



and the evolution of this constraint is subject to the linear 
homogeneous equation 



az„ 



Rfii^Z" 







(3) 



which can be easily obtained from Q allowing for the 
contracted Bianchi identities. Here again, by allowing 
for non-zero values of the extra four- vector Z^, one is 
extending the solution space. But now, as we will see in 
the following section, one is using all the field equations 
to evolve the pair {g^^^Z^}: no equation is taken out 
of the system and, as a result, general covariance is not 
broken. The initial metric g^i, can be taken to be the 
one arising from analytical approximations and the initial 
four-vector Z^ can be taken to vanish without any kind 
of inconsistence: one can even use the evolving values of 
Z^ during the calculation as a good covariant indicator 
of the quality of the approximation. 

Notice that Q plays here the role of the subsidiary 
system. It is adjusted ab initio, without any parameter 
fine-tuning, because light speed is the only characteristic 
speed in Constraint deviations, that is non- vanishing 
values of Z^ will then propagate to the boundaries. The 
fact that our constraints (|2Jl are algebraic will greatly 
simplify the task of providing outgoing boundary con- 
ditions that let the constraint deviations get out of the 
numerical grid [23i | 

Besides these considerations, there are other important 
theoretical issues that we will address in this work. The 
first one is a thorough analysis of the hyperbolicity of the 
evolution system. This is more or less straightforward 
for the first order version of the system, as discussed in 
Appendix B, but it is not so well known in the case of the 
second order version, discussed in Appendix A, where we 
have used the results of Kreiss and Ortiz 26] that recently 
shed light on this issue, which is crucial to discuss the 
well-posedness ,27] of the evolution system (see also 
for similar results for the BSSN system). 

The second theoretical point that we want to stress 
here is that the Z4 system ^ is not just one more 
hyperbolic formalism to be added to the long list. As 
far as it is the only general covariant one, the question 
arises whether the existing non-covariant hyperbolic for- 
malisms |i|l3|ll|lllil|ll|ll|ll|lElIl can be 



recovered from by some "symmetry breaking" mech- 
anism. We have extended in this sense a previous work 
|29| where the deep relationship between the more widely 
used hyperbolic formalisms was pointed out. A partial 
symmetry breaking mechanism is presented in section II 
for recovering the second order systems 0, 0| and for 
the first order systems containing additional dynamical 
fields 9, 10] in section III. These sections are followed 
by another one containing numerical simulations that 
have been proposed recently js^l as standard test-beds 
for Numerical Relativity codes. A more general symme- 
try breaking mechanism is proposed in Appendix C to 
recover first order formalisms which do not contain addi- 
tional dynamical fields El El El El El El . 



II. 3+1 EVOLUTION SYSTEMS 

The general-covariant equations ^ can be written in 
the equivalent 3+1 form 2^ (Z4 evolution system) 

{dt - Cp) 7^, - -2 a (4) 
[dt - Cp) = -V,aj + a ['^^^R,, + V,Z, + W,Z, 



2Kl + {irK - 29) K,, 
Sij + ^{tr S- t) 7y ] 



(5) 
■K 
(6) 



{dt -Cp)e = - + 2 VfeZ*^' + (tiK - 2 9) trK 

- tr{K'^)-2 Z^ak/a-2T] 
{dt - Cf3) Z, ^ a [V, {K,' ~ 6,HtK) + 8,0 

- 2K^ Z,-Qa,/a~ S^-^ (7) 

where we have noted 

9 = a r = 87ra^ r°°, S*, = 87ra T°, = 87r T„ . 

(8) 

In the form H4I7() . it is evident that the Z4 evolution 
system consists only of evolution equations. The only 
constraints lj2Jl, that can be translated into: 
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0, 



z. 



0, 



(9) 



are algebraic so that the full set of field equations 
is actually used during evolution. This is in contrast 
with the ADM evolution system P|, which can be recov- 
ered from H4I7I) by imposing The first two equations 
H4I5|I would transform into the well known ADM evolu- 
tion system, whereas the last two equations (|6l7f) would 
transform into the standard energy and momentum con- 
straints, namely 



^^^R + tr^K -tY{K'^) = 2t 
Vj {K^ - 5,HrK) = 



(10) 
(11) 



In the "free evolution" ADM approach 't'I, both ifTUI) 
and Hll|) were taken out of the evolution system: they 
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were imposed only on the initial data. This was consis- 
tent because (|lUlllf) are first integrals of the ADM evo- 
lution system, but one can not avoid violations of 1)1011111 
due to errors in numerical simulations or approximated 
initial data, as stated before. As a result, numerical sim- 
ulations will deal as well with extended solutions. The 
main difference with the Z4 case, aside from covariance 
considerations, is that in the Z4 case the quantities de- 
scribing constraint deviations are included into the evo- 
lution system. 

One can also ask in this context what happens if, in- 
stead of imposing of the full set one imposes the 
single condition 

9 = (12) 

obtaining a system with only three supplement ary dy- 
namical variables Zi of the kind determined in [29j (Z3 
system): the one corresponding to the parameter choice 

1.1 ^ 2, V ^ n ^ Q (13) 

(we follow the notation of p^). 

One can easily understand two of the three conditions 
((T^ . namely 

fj, = 2, V = n, (14) 



equivalent systems (Z3 evolution systems), namely: 

idt-Cf3)j,, = ~2aK,j (18) 
{dt - Cp)K,j = -V,aj + a f^R.j + V,Zj + VjZ, 

- 2Kl + trK K,^ - S,, + i(tr5 - t)7,,] 

- ^a[^^'^R + 2y-Z + iT^K -iT{K'^) 

-2(a-iafe)Z^ - 2r] 7,, (19) 
{dt - Cp)Z, = a [Vj {K,' - 6^ trK) - 2K,^Z, - 5,] 

(20) 

where we have suppressed the tilde over Kij , allowing for 
the vanishing of Q. 

The resulting system (I18I20|) is quasiequivalent (equiv- 
alent principal parts) to the 'system A' in rcf. |2^ veri- 
fying the 'physical speed' requirement (|14|l . As we have 
already mentioned, it follows that the particular case 

n = I (21) 

is quasiequivalent to the BSSN system 0,^2- The sys- 
tem (I18I21|I can be decomposed into trace and trace-free 
parts 

e^^=7'/', 7., =e-4^7y (22) 



because this amounts to the 'physical speed' requirement 
for the degrees of freedom not related to the gauge 
and nothing else can arise from the general covariant 
equations ^ which are at our starting point. But values 
of the remaining parameter n other than zero would be 
very interesting. In particular, the choice 
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M = 2, jy = n = -. (15) 

would lead to the evolution system which is quasiequiv- 
alent (equivalent principal parts |29|) to the well known 
BSSN system [Hill- 

At this point, let us consider the following recombina- 
tion of the dynamical fields 

K^J = K,, - I e 7y (16) 

so that the Z4 system (I4I7II can be written in a one- 
parameter family of equivalent forms just by replacing 
everywhere 

K,, ^ if,, + |e7„-. (17) 

This kind of transformations leave invariant the solution 
space of the system (it is actually the same system ex- 
pressed in a different set of independent fields). But if 
the suppression of the 8 field H12|) is made after the re- 
placement ()17|l , one gets a one-parameter family of non- 



K = 7^^' K,, , A, = * {K,j - 3 ^ ) (23) 

f . = -7.fc 7 J + 2 Z, (24) 

to follow the correspondence with BSSN more closely. It 
must be pointed out, however, that one does not get in 
this way the original BSSN system: there is actually one 
difference in the lower order terms (only the principal 
parts are equivalent). The difference is in the term of the 
form 

-f^afe Z'= 7y (25) 

in the evolution equation l|19|l . which is missing in the 
orig inal BSSN system This lower order term is 

needed for consistency with the general covariant equa- 
tions 

We have seen then how the widely used ADM and 
BSSN systems can be obtained from the more general 
Z4 formalism. The equivalence transformation (|16|l plays 
the crucial role because suppressing the Q field H12|l pro- 
duces a sort of symmetry breaking: different values of 
the parameter n will lead to evolution systems that can 
no longer be transformed one into another once the set of 
dynamical fields is reduced by the disappearance of O. It 
can be regarded as a partial symmetry breaking mecha- 
nism for the original equations (|4I6|I . The terms "partial" 
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refers to the fact that only the quantity Q is suppressed, 
while the are kept into the system (|18I20|I . A complete 
symmetry breaking mechanism is discussed in Appendix 
C. 

In section IV, we present the results of some test-bed 
simulations for the ADM and Z4 systems. We have con- 
sidered for simplicity only vacuum space-times with the 
time coordinate conditions 

{dt - Cfi) In a - a [fivK ~ AG] (26) 

which are a further generalization of the one proposed in 
where 

/ = 1, A = 2. (27) 

This two-parameter family of coordinate conditions is 
very interesting from the point of view of Numerical Rel- 
ativity applications. But it is also interesting from the 
theoretical point of view, because it provides the oppor- 
tunity to apply the recent results of Kreiss and Ortiz 
on the hyperbolicity of the ADM system in a wider con- 
text. In Appendix A, we will use the same formulation 
(see ref. j23| for more details) to study the hyperbolic- 
ity of the Z4 system with the two-parameter family of 
dynamical gauge conditions H26|l . 

III. FIRST ORDER SYSTEMS 

A first order version of the Z4 evolution system 14I7|) 
can be obtained in the standard way by considering the 
first space derivatives 

Ak = ak/a, Dkij = ^ dki.j (28) 

as independent dynamical quantities with evolution 
equations given by 

dtAk + dk[ a{f tvK-Xe)] = (29) 



^tDk^J + dk[aK,j] = (30) 

(we will consider in what follows the vanishing shift case 
for simplicity), so that the full set of dynamical fields can 
be given by 

u = {a, 7y, Kij, Ak, Dktj, 6, Zk} (31) 

(38 independent fields). 

Care must be taken when expressing the Ricci tensor 
in lO in terms of the derivatives of Dkij , because 
as far as the constraints (|28|l are no longer enforced, the 
identity 

Qj-Dgij — dgDrij (^^) 

can not be taken for granted in first order systems. As a 
consequence of this ordering ambiguity of second deriva- 
tives, the principal part of the evolution equation can 



be written in a one-parameter family of non-equivalent 
ways, namely 

dtK.j + dk [a X%] = ... (33) 

+ 1 (A, + D,/ - 2Z,) + i S';iA, + A/ - 2Z,) 

(34) 

so that the parameter choice C = +1 corresponds to the 
standard Ricci decomposition 

— dk r'^ij — di T'^kj + ^^rk^'^ij — ^'^ri^^kj (35) 

whereas the opposite choice ( = —I corresponds to the 
de Donder-Fock [3lll3^ decomposition 

'^^^Rij — —dk D^ij + d(i Tj-fk^ — 2D/^Dkij 

+ '^D^^iDrsj — ^irs^f^ — ^rij^^^k (36) 

which is most commonly used in Numerical Relativity 
codes. Note that this ambiguity does not affect to the 
principal part of eq. ©, namely 

dtQ + dk[aV''] = ... (37) 

where we have noted 

Vk EE Dk/ - D\k - Zk. (38) 

We are now in position to discuss the hyperbolicity of 
this first order version of the Z4 systems. This is done in 
a straightforward way in Appendix A. 

In order to compare the new first order system with 
the Bona-Masso ones 0, ^| we could either apply here 
again the recombination H16|l followed by the suppression 
(|12|l of the 8 field or we could take directly a first order 
version of the Z3 system (I18I2U|) . Eqs. H18I2UI30|I do not 
change, but eqs. I|29I33I34|) are replaced in this case by 

dtAk + dk [a / tiK] = (39) 
dtK.j + dk [a A^,] ^ ... (40) 

^ -r%-|vH, (41) 

+ ^SUA, + D,/ - 2Z,) + ^S^{A, + A/ - 2Z,) 

The full Bona-Masso family of evolution equations is re- 
covered for the ^ = — 1 case, where ()41|l can be written 
as 

Af^. ^ D% - (42) 
+ l<5f (A, - A/' + 2y,) + ls^{A, - A."- + 2Vi) 
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with Vk defined by 

In the following section, we will compare the behavior 
of both families in numerical simulations. To this end, 
we will also consider the first order version of the ADM 
system, which can be obtained from the previous ones 
just by suppressing the Zi eigenfields. 

Z, ^ 0. (43) 

Before proceeding to the test section, let us just men- 
tion, that the same game of recombining the field with 
Kij (I17|l before suppressing it can also be played with 
the Zk fields and Dkij in first order systems. As stated 
before, this will provide a complete symmetry breaking 
mechanism. We will do that in Appendix C, where we 
will show how the well known KST system can also 
be recovered in that way from the Z4 framework dis- 
cussed in this paper. 



IV. TESTING SECOND AND FIRST ORDER 
SYSTEMS 

We will present in this section a couple of numerical ex- 
periments which have been suggested very recently 30] as 
standard test-beds for Numerical Relativity codes. Our 
philosophy is that all the tests could be done "out of the 
box" , using well known numerical methods and the equa- 
tions that are fully presented here: anyone should be able 
to reproduce our results without recourse to additional 
information. 

We will use the standard method of lines as a fi- 
nite differencing algorithm, so that space and time dis- 
cretization will be dealt separately. Space differencing 
will consist of taking centered discretizations of deriva- 
tives in our 3D grid. We use the standard centered stencil 
for first derivatives and we make sure that second deriva- 
tives, when needed, are coded also as centered deriva- 
tives of these first derivatives, even if it takes up to five 
point along every axis. In order to avoid boundary effects, 
the grid has the topology of a three-torus, with periodic 
boundaries along every axis. The time evolution will be 
dealt with a third order Runge-Kutta algorithm. The 
time step dt is kept small enough to avoid an excess of 
numerical dissipation that could distort our results in 
long runs. 



A. Robust stability test 

Let us consider a small perturbation of Minkowski 
space-time which is generated by taking random initial 
data for every dynamical field in the system. The level 
of the random noise must be small enough to make sure 
that we will keep in the linear regime even for a thousand 
of crossing times (the time that a light ray will take to 
cross the longest way along the numerical domain) . This 
is in keeping with the theoretical framework of Appendix 
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FIG. 1: The maximum of (the absolute value of) trK is 
plotted against the number of crossing times in a logarithmic 
scale. The initial level of random noise remains constant dur- 
ing the evolution in the case of strongly hyperbolic systems 
(only the second order Z4 system is shown here for clarity). 
In the case of weakly hyperbolic systems, like the ADM sec- 
ond order system ADM-2 or its first order version ADM-l, 
a linear growth is detected up to the point where the codes 
crash. Notice that the second order version is more robust, 
an order of magnitude, than the first one. The simulations 
are made with 50 grid points with dt = 0.03 da;. 



A, where only linear perturbations around the Minkowski 
metric are considered. 

We have plotted in Fig. our results for the standard 
harmonic case (I27II . We see the expected polynomial (lin- 
ear in this case) growth [23| of the weakly hyperbolic 
ADM system. Notice that modifications of the lower 
order terms (the ones not contributing to the principal 
part) could lead to catastrophic exponential growth, re- 
vealing an ill-posed evolution system 27]. In this paper, 
however, we will limit ourselves to discussing the linear 
regime as an hyperbolicity test for the principal part of 
the system. In this sense, the linear growth of the ADM 
plots in Fig. ^ confirms the weakly hyperbolic character 
of the ADM system. 

The Z4 system shows instead the no-growth behavior, 
independent of the time resolution (see Figs.^lU, which 
one would expect from a strongly hyperbolic system. The 
same qualitative behavior is shown by the corresponding 
Z3 systems in (|18I20|I . including the one with n = 4/3 
that is quasiequivalent to the BSSN system. 

We also show in Fig.^jthe same results, but distorted 
by using too much numerical dissipation: the time evo- 
lution here is dealt with the second order ICN method 
rather than the third order Runge-Kutta of Fig. Af- 
ter hundreds of crossing times, the numerical dissipation 
manages to curve the linear growing of ADM and the 
noise level goes down in the Z4 case. This is just a nu- 
merical artifact, because in the linear regime there is no 
physical damping mechanism for strongly hyperbolic sys- 
tems in a three-torus, where periodic boundary condi- 
tions do not allow propagation outside the domain. This 
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time time 



FIG. 2: Same as Fig. but with less time resolution 
{dt — 0.06 da; with the same da;). There is a slight amount 
of dissipation that delays the crashing of the ADM codes; 
this is especially visible for the second order version ADM-2, 
which keeps being more robust than its first order counterpart 
ADM-1. The behaviour of the Z4 code keeps unaffected. 



is why we will use here third order Runge-Kutta instead 
of the ICN method proposed in [3(tI |. 

In Fig. O we explore parameter space in the (/,A) 
plane. If we interpret the constant behavior in Fig. ^ 
as revealing a strongly hyperbolic system and the poly- 
nomial growth (linear in this case) in Fig. ^ as revealing 
a weakly hyperbolic system, the results of our numerical 
experiment fully agree with the theoretical results pre- 



FIG. 3: Array of results of numerical experiments in the 
gauge parameter plane (/,A), by using the Z4 system. A 
triangle stands for linear growth of noise (weak hyperbolic- 
ity), whereas a cross stands for a constant noise level (strong 
hyperbolicity). This is consistent with the strong hyperbolic- 
ity requirements predicted in Appendix B: either f — 1 and 
A = 2, or / / 1 and / > 0. 

4 X A X X X 

3 X A X X X 

2 A X X X X X 

1 A X A X X X 

" X A X X X 

-1 A X A X X X 



2 A X A- 



FIG. 4: Same as Fig. El but using the second order ICN 
method to evolve in time instead of a third order Runge-Kutta 
algorithm. Numerical dissipation is severely distorting the 
plots, by masking the linear growth in the weakly hyperbolic 
case and dramatically reducing the initial noise level in the 
strongly hyperbolic case. Notice than both dt and da; are 
the same as in Fig. Q and we are using also the same space 
discretization algorithm: only the time evolution method has 
changed. 



sented in Appendix A. 



B. Gowdy waves 

In order to test the strong field regime, let us consider 
now the Gowdy solution |2J|, which describes a space- 
time containing plane polarized gravitational waves (see 
also for an excellent review of these space-times as 
cosmological models). The line element can be written 
as 

ds^ = e«/2 {-dt^+dz^)+t (e^ dx^+e-^ dy^) (44) 

where the quantities Q and P are functions of t and z 
only and periodic in z, so that H44|l is well suited for 
finite difference numerical grids with periodic boundary 
conditions along every axis. Following , we will choose 
the particular case 

P = Jo(27ri) cos(27rz) (45) 
Q = 7rJo(27r)Ji(27r) 27r<Jo(27ri)Ji(27rt)cos^(27rz) 
-I- 2TTh^[Jo^{2TTt) + Ji^{2TTt) - Jo^(27r) - Ji^(27r)] 

(46) 



so that it is clear that the lapse function 



(47) 



X > 
4 / 



is constant e very where at any time at which Jo(27rio) 
vanishes. In [23 the initial slice t = to was chosen for the 
simulation of the collapse, where 27rto is the 20-th root 
of the Bessel function Jo, i-e. to — 9.88. 
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FIG. 5: Time evolution of (the maximum value of) the lapse 
function a in a collapsing Gowdy space-time (harmonic slic- 
ing). Notice that the harmonic time coordinate t is not the 
proper time and it does not coincide with the number of cross- 
ing times, due to the collapse of the lapse, which is visible 
here, by a 10~^ factor. 



FIG. 7: Convergence test for the Z4 code. The three lines 
correspond to 50, 100 and 200 grid points. The quantity O 
itself is a direct measure of the error. In that logarithmic scale 
differences of log 4 correspond to dividing by four the error 
when doubling the resolution. This second order convergence 
rate is clearly shown in the figure. 



Let us now perform the following time coordinate 
transformation 



t = to e 



To 



y^3/4„Q(to)/4 



472 



(48) 



so that the expanding line element (|44|l is seen in the 
new time coordinate t as collapsing towards the t = 
singularity, which is approached only in the limit r — > oo . 
This "singularity avoidance" property of the r coordinate 
is not surprising if one realizes that the resulting slicing 
by r = constant surfaces is harmonic [s^ . 

This means that we can launch our simulations starting 
with a constant lapse ag = 1 at r = (t = to) with the 



FIG. 6: The quantity B is plotted as an indicator of the 
accumulated error of the simulations for the ADM, Z4 and 
Z3-BSSN second order forms. Even in this logarithmic scale, 
it can be clearly seen how the Z4 and Z3-BSSN codes perform 
much better than the ADM one: error differs by one order of 
magnitude at r ~ 1000. The Z3-BSSN code gets closer to the 
Z4 one in the oscillatory phase (up to r ~ 2000). 



ADM , 




gauge parameter choice / = 1 (which means also A = 2 in 
the Z4 case). Notice that the harmonic time coordinate 
T is not the proper time and it does not coincide with the 
number of crossing times, due to the collapse of the lapse. 
Remember also that the local value of light speed (proper 
distance over coordinate time) is ±a. Even though in our 
plots T goes up to 10000, the hght ray manages to cross 
the domain in the z-direction only to ~ 9.88 times, as it 
follows from the original form (|44|l of the line element. 

We plot in Fig. El the maximum values of the lapse 
function as time goes on, measured in terms of the har- 
monic time coordinate r. Notice the huge magnitude of 
the dynamical space we are covering, as a goes down (col- 



FIG. 8: Same as Fig. El but with the first order versions of 
both the ADM and Z4 codes, which behave in the same way 
as their second order counterparts. The Z3-BM code gets here 
closer the ADM one in the oscillatory phase (up to r ~ 2000), 
in contrast with the behavior of the Z3-BSSN code in Fig.jSl 



ADM , 



Z3-BM 



1000 10000 




10000 



8 



lapse of the lapse) by the factor of one billion during the 
simulation. This is a real challenge for numerical codes 
and all of them are doing quite well until r = 1000. The 
behavior at later times is dominated by the lower order 
terms: coordinate light speed (±a) is so small that the 
dynamics of the principal part is frozen and care must 
be taken to avoid too big time steps. We have not seen 
any of the codes crashing even in very long simulations 
(up to T = 60.000) when the size of the time step is kept 
under control. 

In Fig.ini we use the quantity 0, as computed from 
to monitor the quality of the simulation both in the Z4 
case (|4I7|I , where is a dynamical field, and in the other 
two second order cases (ADM and Z3-BSSN), where O 
is no longer a dynamical quantity but can still be used 
as a good measure of the error in the simulation. This 
makes easier to perform convergence tests (second order 
convergence is shown in Fig. [71). Notice that our Z3-BSSN 
code performs here much better than the original BSSN 
one J as reported in j^Ol ■ This is mainly due to the fact 
that we are not using here the conformal decomposition 
H22I24|I which is at odds with the structure of the line 
element (glj. This is why we talk here about Z3-BSSN 
()18li>l|) instead of simply BSSN. 

The same kind of comparison is made in Fig.|Slfor the 
first order version of the ADM and Z4 codes, which show 
the same behavior than their second order counterparts 
in Fig. El The third plot corresponds to the Z3 version of 
the Bona-Masso code, which can be obtained from H39t 
I41|l . with the parameter choice — —1, n — +1. Notice 
that in the oscillatory phase (up to r ~ 2000) the Z3-BM 
code tends to behave like the ADM one, whereas in Fig. El 
the Z3-BSSN code tends to behave more like the Z4 one. 



FIG. 9: Same as Fig. |S1 but with different values of the 
parameter n arising from the symmetry breaking mechanism 
in the Z3 codes. The differences show up in the collapse final 
phase (starting at r ~ 2000). Notice that the value n = 4/3 
corresponds to the Z3-BSSN case in Fig. |^ whereas the case 
n = 1 corresponds to (the second order version of) the Z3-BM 
case in Fig. |H1 First order versions (not shown) behave in the 
same way as their second order counterparts shown here 



10+1" n 




Notice that the Z3 versions, when compared with both 
ADM and Z4, show a different behavior after the oscilla- 
tory phase: grows at a much lower rate or even starts 
going down. This behavior is due to the extra terms that 
appear in the evolution equations for Kij after the sym- 
metry breaking, which can be controlled by the parame- 
ter n introduced in (|17|l . This point is clearly shown in 
Fig. [51 where different choices of n produce different be- 
havior in the final collapse phase (starting at r ~ 2000). 
This shows the relevance of the recombination between 
the dynamical fields in numerical applications, as pointed 
out in jlj| . Further details and other numerical tests can 
be found in the webpage at http://stat.uib.es. 



Appendix A: Hyperbolicity of the Second Order Z4 
System 

Let us consider the linearized version of the Z4 system 
(I4I7II around Minkowski space-time in order to study the 
propagation of a plane wave in that background: 



7y = 5,,+2e'^-'' %{uj,t) (A.l) 

a = iH-e'"-^ d(a;,i) (A.2) 

= i uj e"^-"" ky{uj,t) (A.3) 

e = iuje"^-"" Q{uj,t) (A.4) 

Zk = ioje"^-"" Zk{uj,t) (A.5) 

where we will take for simplicity /?* = and 

ujk ^ uj Uk , S"-^ Hi Uj = 1 (A. 6) 

The Z4 system reads then 

9t7y = -iui Kij (A. 7) 

dta = w [/ trk - A 8] (A.8) 

dtQ ^ ^iio [tr7 - 7"" - Z"] (A.9) 
dtZu = -iio K [trk - 6) - k-^] (A. 10) 

dtkij = -iio Xij (A. 11) 



where we have noted 

\j = % + n, rij (d + tr^) 

- (7; + Z./) - n, (77 + Z,) (A.12) 

and where the symbol n replacing an index means the 
contraction with n^. It can be also expressed in matrix 
form, namely 

u = {a,^ij,kij,Q,Zk) (A. 13) 

dtu = -iuj Ail (A. 14) 



The spectral analysis of the characteristic matrix A 
provides the following list of eigenvalues and eigenfields 
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a) Standing eigenfields (zero characteristic speed) 

a - /tr7 + A(tr7 - 7"" - Z") , T^ + Z^ (A.15) 

where the symbol _L replacing an index means the 
projection orthogonal to n^. 

b) Light-cone eigenfields (characteristic speed ±1) 



e ± [tr7 - f 



ncnts of a tensor: 



[kij - [trk - 29) n., nj\ 



± [Ay 

c) Gauge eigenfields (characteristic speed ±\ff) 
± 



/-I 



(tr7 - 7"" - Z") 



(A.20) 



2/-A 
/-I 



From HA.15IA.20p we can easily conclude |23 

i) All the characteristic speeds are real (weak hyper- 
bolicity at least) if and only if / > 0. 

ii) In the case / = 0, the two components of the gauge 
pair (|A.20p are not independent, so that the total 
number of independent eigenfields is 16 instead of 
17 required for strong hyperbolicity. 

iii) The case / = 1 (harmonic case) is special: 



• if A 7^ 2, then the gauge pair (|A.20p . which 
can be previously rescaled by a (/ — 1) factor, 
is equivalent to IjA.lSp . so that one has only 
15 independent eigenfields 

• if A = 2, then the quotient can take any 
value, reflecting the degeneracy of the gauge 
and light cone eigenfields. One can then re- 
cover the full set of 17 independent eigenfields 
(strong hyperbolicity). 

iv) In all the remaining cases (/ > 0, / 7^ 1), the 
system is strongly hyperbolic, as we can recover 
the full set of 17 independent eigenfields. 



Appendix B: Hyperbolicity of the First Order Z4 
System 

The principal part of the first order Z4 evolution sys- 
tem (|4I7I26I29I3(J|I can be written as (vanishing shift 



dtlij = ■ • ■ , dta^ ... 

dtO + dk[aV'']^... 

dtZ, + dk [a (S^itvK - 9) - ^;)] = 

dtAk + dk [a {ftvK - A9)] = . . . 

dtDkij + dk [a K^j] = . . . 



(A.16) 






(A.17) 


where 




(A.18) 


Ay - 




compo- 


+ 
+ 




(A.19) 


Vk = 





(B.l) 
(B.2) 
(B.3) 
(B.4) 
(B.5) 
(B.6) 



2K~ 



^rk 



(B.7) 
(B.8) 



Now, if we consider the propagation of perturbations 
with wavefront surfaces given by the unit (normal) vector 
rii, we can express (|B.llB.6ll in matrix form 

a~^dt u + A(m) v!' dkU = ... , 



where 

u = {a, 7y-, Kij, Ak, Dkij, 9, Zk} 



(B.9) 
(B.IO) 



(notice that derivatives tangent to the wavefront surface 
play no role here). 

A straightforward analysis of the characteristic matrix 
A(m) provides the following list of eigenfields: 

a) Standing eigenfields (zero characteristic speed) 

a, 7y, A^, D^,j, Ak - fDk + XVk (B.ll) 

(24 independent fields), where the symbol _L replac- 
ing an index means the projection orthogonal to nf. 



D 



Dkij - nkffDr 



(B.12) 



b) Light-cone eigenfields (characteristic speed ±1) 



L ij = [Kij - niJij ivK] 
± [A"y - niUj tr A'^ 



(B.13) 
(B.14) 



(12 independent fields), where the symbol n replac- 
ing the index means the contraction with Ui 



Ay- = nk Xfj. 



(B.15) 



c) Gauge eigenfields (characteristic speed ±v7) 

2 - A 



tr K 



± 



./-I 

2/ -A 



9 



/-I 



(B.16) 
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From ljB.lllB.16ll we can easily conclude that 

i) All the characteristic speeds are real (weak hyper- 
bolicity at least) if and only if / > 0. 

ii) In the case f — 0, the two components of the pair 
(jB.16|l are not independent, so that the total num- 
ber of independent eigenfields is 37 instead of 38 
required for strong hyperbolicity. 

iii) The case f=l (harmonic case) is special: 

• if A 7^ 2, then the pair of fields ljB.16|) is the 
same as (|B.14ll . so that one has only 36 inde- 
pendent eigenfields 

• if A = 2, then the quotient can take any 
value due to the degeneracy of the gauge and 
light eigenfields. One can then recover the full 
set of 38 independent eigenfields (strong hy- 
perbolicity) . 

iv) The first order Z4 system described by ljB.llB.6|l is 
strongly hyperbolic in all the remaining cases (/ > 
0, /^l). 

Notice also that the special (harmonic) case / = 1, 
A = 2 has been shown in |25| to be symmetric hyperbolic 
for the parameter choice ^ = — 1. The corresponding 
energy function can be written as 

E = K'^K.j + X'^'^Xmj + {trK - 26)^ + A'^Ak 
+ (A'' - D''% + 2V''){Ak- Dk/ + 2Vk) (B.17) 

but notice that this expression is far from being unique. 
For instance, allowing for ljB.ll|) . the last term in ljB.17|) 
could appear with any arbitrary factor. 

Appendix C: Recovering the KST systems 

Let us start with the first order Z4 evolution system 
(|4l7l26l29l30f) where the principal part is given by IjB.lt 
IB.8|I . Now let us follow the two step 'symmetry breaking' 
process, namely 

i) Recombine the dynamical fields Kij^ D^ij with O 
and Zi in a linear way, 

K., - i^., - 2 , (C.l) 

dkij = 2Dkij + f] lk(iZj) + X Zkli] , (C.2) 

where we have used the notation of Ref. 0|, re- 
placing only their parameter 7 by ~n/2 for con- 
sistency. Notice that ljC.llC.2|) is generic in the 
sense that it is the most general linear combination 
that preserves the tensor character of the dynami- 
cal fields under linear coordinate transformations. 



ii) Suppress both 6 and Zi as dynamical fields, namely 

9 = 0, Z, = . (C.3) 

In that way, the principal part (|B.1IB.8|I becomes 

dt%j = ... , dta^... (C.4) 
dtAk + dk [a firk] = (C.5) 

dtdk:, + dr[a{2 8lk,,-x{kl-5ltvk)^,, 

+ iniK:{k^,)~5^,)^^K)}]^--- (C.6) 

dtk,, + dk [a \%] = . . . (C.7) 

Z A ij — LI ij ^ \^kr Jiij 

+ (A, + id„'-) + Sf (A, + id,/) (C.8) 

for the reduced set of variables 

u = {a, 7ij, ki-j, Ak, dkij}. (C.9) 

This provides a "dynamical lapse" version of the 
KST evolution systems. In order to recover the original 
"densitized lapse" version, one must in addition integrate 
explicitly the dynamical relationship H26|l between the 
lapse and the volume element (remember that now O = 
0). It can be easily done in the case 

/ = 2 (T = constant, (C.IO) 

namely 

5t(a7-") = 0, (C.ll) 

so that the value of a can be defined in terms of 7 for 
every initial condition. The same thing can be done with 
Ai and di, so that 

A, = ad,/ + ... (C.12) 
and the set of dynamical fields is then further reduced to 

u = {jij, Ktj, dkij}. (C.13) 

The principal part of the evolution system is then given 
by (we suppress the tildes over the Kij) 

dti^j = ... (C.14) 

dtdktj + dr[a{2 8l K,j - x {K/ - 51 trK)-^.^ 

+ V lk(^{K^f) - 5]^^^K) }] = ... (C.15) 

dtK,, + dk [a Af,] = . . . (C.16) 
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+ ^^^{S^d,/ + S^d,/) (C.17) 

which corresponds precisely to (the principal part of) 
the original KST system . 
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